set more off
use Samuels_Thomson_PoP_Data.dta, clear
set scheme s2mono
xtset cown year

************************
* ANALYSIS STARTS HERE *
************************


*****************
* Tractors only *
*****************
* Polity
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha, fe vce(cluster cown)
est store tractors_polity

margins, at(L4.z2_l_ln_tractors_ha=(-1(0.25)1))
marginsplot, recast(line) recastci(rline) plot1opts(lwidth(medthick)) ciopts(lwidth(medium) lpattern(longdash)) ///
addplot(kdensity z2_l_ln_tractors_ha if e(sample) & z2_l_ln_tractors_ha>-1  & z2_l_ln_tractors_ha<1, lpattern(dot) lwidth(medthick) yaxis(2)) ///
title("Model 1.1", pos(12)) ///
xtitle("Tractors") xlabel(-1(0.25)1) ///
ytitle("Polity", axis(1)) ///
ysize(5) xsize(6) graphregion(color(white)) bgcolor(white) ///
legend(off)
graph save polity_tractors, replace

* Unstandardized to aid interpretation
xtreg polity_s_2010 L4.z2_l_ln_tractors_ha, fe vce(cluster cown)
margins, at(L4.z2_l_ln_tractors_ha=(-1(0.25)1))

* VDem Polyarchy
xtreg z2_v2x_polyarchy L4.z2_l_ln_tractors_ha, fe vce(cluster cown)
est store tractors_polyarchy
* VDem Liberal Democracy
xtreg z2_v2x_libdem L4.z2_l_ln_tractors_ha, fe vce(cluster cown)
est store tractors_libdem


*****************
* Plus Controls *
*****************
* Polity
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store tractors_polity_controls
estat ic

gen model2_sample = 0
replace model2_sample = 1 if e(sample)

margins, at(L4.z2_l_ln_tractors_ha=(-1(0.25)1))
marginsplot, recast(line) recastci(rline) plot1opts(lwidth(medthick)) ciopts(lwidth(medium) lpattern(longdash)) ///
addplot(kdensity z2_l_ln_tractors_ha if e(sample) & z2_l_ln_tractors_ha>-1  & z2_l_ln_tractors_ha<1, lpattern(dot) lwidth(medthick) yaxis(2)) ///
title("Model 1.2", pos(12)) ///
xtitle("Tractors") xlabel(-1(0.25)1) ///
ytitle("Polity", axis(1)) ///
ysize(5) xsize(6) graphregion(color(white)) bgcolor(white) ///
legend(off)
graph save polity_tractors_controls, replace

* Compare to Model w/o Tractors
xtreg z2_polity_s_2010 ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc if e(sample), fe vce(cluster cown)
est store polity_controls
estat ic

* LR Test of Models 
lrtest polity_controls tractors_polity_controls, force stats


* Non-standardized to aid substantive interpretation
xtreg polity_s_2010 L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2, fe vce(cluster cown)
margins, at(L4.z2_l_ln_tractors_ha=(-1(0.25)1))

****************
* VDem Polyarchy 
xtreg z2_v2x_polyarchy L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store tractors_polyarchy_controls
estat ic


margins, at(L4.z2_l_ln_tractors_ha=(-1(0.25)1))
marginsplot, recast(line) recastci(rline) plot1opts(lwidth(medthick)) ciopts(lwidth(medium) lpattern(longdash)) ///
addplot(kdensity z2_l_ln_tractors_ha if e(sample) & z2_l_ln_tractors_ha>-1  & z2_l_ln_tractors_ha<1, lpattern(dot) lwidth(medthick) yaxis(2)) ///
title("Model 1.3", pos(12)) ///
xtitle("Tractors") xlabel(-1(0.25)1) ///
ytitle("VDem Polyarchy", axis(1)) ///
ysize(5) xsize(6) graphregion(color(white)) bgcolor(white) ///
legend(off)
graph save polyarchy_tractors_controls, replace

* Compare to Model w/o Tractors (VDem Polyarchy)
xtreg z2_v2x_polyarchy ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc if e(sample), fe vce(cluster cown)
estat ic


************************
* VDem Liberal Democracy
xtreg z2_v2x_libdem L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store tractors_libdem_controls
estat ic

margins, at(L4.z2_l_ln_tractors_ha=(-1(0.25)1))
marginsplot, recast(line) recastci(rline) plot1opts(lwidth(medthick)) ciopts(lwidth(medium) lpattern(longdash)) ///
addplot(kdensity z2_l_ln_tractors_ha if e(sample) & z2_l_ln_tractors_ha>-1  & z2_l_ln_tractors_ha<1, lpattern(dot) lwidth(medthick) yaxis(2)) ///
title("Model 1.4", pos(12)) ///
xtitle("Tractors") xlabel(-1(0.25)1) ///
ytitle("VDem Lib Dem", axis(1)) ///
ysize(5) xsize(6) graphregion(color(white)) bgcolor(white) ///
legend(off)
graph save libdem_tractors_controls, replace

* Compare to Model w/o Tractors (VDem Lib Dem)
xtreg z2_v2x_libdem ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc if e(sample), fe vce(cluster cown)
estat ic


**********************
* Year Fixed Effects *
**********************
* Polity
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc i.year, fe vce(cluster cown)
est store tractors_polity_years

* VDem Polyarchy
xtreg z2_v2x_polyarchy L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc i.year, fe vce(cluster cown)
est store tractors_polyarchy_years

* VDem Liberal Democracy
xtreg z2_v2x_libdem L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc i.year, fe vce(cluster cown)
est store tractors_libdem_years


*********************************
* Labor-Dependent Ag (Albertus) *
*********************************
* Polity
xtreg z2_polity_s_2010 z2_l_repressive_labor_pc L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store lrag_polity

* Compare magnitude of LR Ag and Mechanization
margins, at(L4.z2_l_ln_tractors_ha=(-1(0.25)1))
marginsplot, recast(line) recastci(rline) plot1opts(lwidth(medium)) ciopts(lwidth(thin) lpattern(longdash)) ///
title("Model 1.6", pos(12)) ///
xtitle("Tractors") xlabel(-1(0.25)1) ///
ytitle("Polity") ///
addplot(kdensity z2_l_ln_tractors_ha if e(sample) & z2_l_ln_tractors_ha>-1 & z2_l_ln_tractors_ha<1, lpattern(dot) lwidth(medthick) yaxis(2)) ///
legend(off) xsize(4) ysize(3) ///
graphregion(color(white)) bgcolor(white)
graph save polity_tractors_3, replace

margins, at(z2_l_repressive_labor_pc=(-1(0.25)1))
marginsplot, recast(line) recastci(rline) plot1opts(lwidth(medthick)) ciopts(lwidth(medium) lpattern(longdash)) ///
addplot(kdensity z2_l_repressive_labor_pc if e(sample) & z2_l_repressive_labor_pc>-1  & z2_l_repressive_labor_pc<1, lpattern(dot) lwidth(medthick) yaxis(2)) ///
title("Model 1.6", pos(12)) ///
xtitle("Labor-Dependent Ag") xlabel(-1(0.25)1) ///
ytitle("Polity", axis(1)) ///
ysize(5) xsize(6) graphregion(color(white)) bgcolor(white) ///
legend(off)
graph save polity_lrag, replace

* Albertus' interactions with civil war, rural unrest
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc l_civilwarprio repressive_labor_civilwarp ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store lrag_civilwar_polity

xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc l10_minorland_natl_avg_excl rl_l10_minorland_natl_avg_excl ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store lrag_reform_polity

* Polyarchy
xtreg z2_v2x_polyarchy z2_l_repressive_labor_pc L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store lrag_polyarchy

* Liberal Democracy
xtreg v2x_libdem z2_l_repressive_labor_pc L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store lrag_libdem


****************************
* Family Farms (Van Hanen) *
****************************
xtreg z2_polity_s_2010 z2_l_famfarms L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store famfarms_polity

margins, at(z2_l_famfarms=(-1(0.25)1))
marginsplot, recast(line) recastci(rline) plot1opts(lwidth(medthick)) ciopts(lwidth(medium) lpattern(longdash)) ///
addplot(kdensity z2_l_famfarms if e(sample) & z2_l_famfarms>-1  & z2_l_famfarms<1, lpattern(dot) lwidth(medthick) yaxis(2)) ///
title("Model 1.7", pos(12)) ///
xtitle("Family Farms") xlabel(-1(0.25)1) ///
ytitle("Polity", axis(1)) ///
ysize(5) xsize(6) graphregion(color(white)) bgcolor(white) ///
legend(off)
graph save polity_famfarms, replace

* Polyarchy
xtreg z2_v2x_polyarchy z2_l_famfarms L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store famfarms_polyarchy

* Liberal Democracy
xtreg v2x_libdem z2_l_famfarms L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store famfarms_libdem



******************
* Other controls *
******************
* Inequality
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha l_siddgini_interp ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)

* Oil Income
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha l_total_oil_income_pc_thou ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)

* Capital Acct Openness
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha l_kaopen ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)



************
* FIGURE 4 *
************
graph combine polity_tractors.gph polity_tractors_controls.gph polyarchy_tractors_controls.gph libdem_tractors_controls.gph ///
polity_lrag.gph polity_famfarms.gph, cols(2) ycommon graphregion(color(white)) xsize(7) ysize(8)


***********
* TABLE 1 *
***********
esttab tractors_polity ///
tractors_polity_controls tractors_polity_years tractors_polyarchy_controls tractors_libdem_controls ///
lrag_polity ///
famfarms_polity, ///
keep(L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc ///
z2_l_famfarms ///
z2_l_log_gdp_per_cap_haber_men_2 ///
z2_l_region_dems_pc) ///
order(L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc ///
z2_l_famfarms ///
z2_l_log_gdp_per_cap_haber_men_2 ///
z2_l_region_dems_pc) ///
coeflabels(L4.z2_l_ln_tractors_ha "Tractors" ///
z2_l_repressive_labor_pc "LD Ag" ///
z2_l_famfarms "Fam Farms" ///
z2_l_log_gdp_per_cap_haber_men_2 "GDP/Cap" ///
z2_l_region_dems_pc "Reg Dem (%)") ///
b(2) scalar(N_clust) aic se nogaps compress ///
star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("Polity" "Polity" "Polyarchy" "Lib Dem" "Polity" "Polity" "Polity" "Polity") ///
addnotes("All models with country fixed effects, clustered standard errors" ///
"Model 5 with year fixed effects")

***********
* TABLE 2 *
***********
esttab tractors_polyarchy tractors_libdem ///
tractors_polyarchy_years tractors_libdem_years ///
lrag_polyarchy lrag_libdem ///
lrag_civilwar_polity lrag_reform_polity ///
famfarms_polyarchy famfarms_libdem, ///
keep(L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc ///
l_civilwarprio repressive_labor_civilwarp ///
l10_minorland_natl_avg_excl rl_l10_minorland_natl_avg_excl ///
z2_l_famfarms ///
z2_l_log_gdp_per_cap_haber_men_2 ///
z2_l_region_dems_pc) ///
order(L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc ///
l_civilwarprio repressive_labor_civilwarp ///
l10_minorland_natl_avg_excl rl_l10_minorland_natl_avg_excl ///
z2_l_famfarms ///
z2_l_log_gdp_per_cap_haber_men_2 ///
z2_l_region_dems_pc) ///
coeflabels(L4.z2_l_ln_tractors_ha "Tractors" ///
z2_l_repressive_labor_pc "LD Ag" ///
l_civilwarprio "Civil Conflict" repressive_labor_civilwarp "LD Ag*Conflict" ///
l10_minorland_natl_avg_excl "Land Reform" rl_l10_minorland_natl_avg_excl "LD Ag*Land Ref" ///
z2_l_famfarms "Fam Farms" ///
z2_l_log_gdp_per_cap_haber_men_2 "GDP/Cap" ///
z2_l_region_dems_pc "Reg Dem (%)") ///
b(2) scalar(N_clust) aic se nogaps compress ///
star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("Polyarchy" "Lib Dem" "Polyarchy" "Lib Dem" "Polyarchy" "Lib Dem" "Polity" "Polity" "Polyarchy" "Lib Dem") ///
addnotes("All models with country fixed effects, clustered standard errors" ///
"Models 2.3-2.4 with year fixed effects")


*******************************
* Additional Robustness Tests *
*******************************

**************
* Time Trend *
**************
gen time = year-1935

* Polity
xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc c.time##c.time##c.time, fe vce(cluster cown)
est store polity_cubic_time

* VDem Polyarchy
xtreg z2_v2x_polyarchy L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc c.time##c.time##c.time, fe vce(cluster cown)
est store polyarchy_cubic_time

* VDem Liberal Democracy
xtreg z2_v2x_libdem L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc c.time##c.time##c.time, fe vce(cluster cown)
est store libdem_cubic_time


********************************
* Use same sample as Model 1.6 *
********************************
xtreg z2_polity_s_2010 z2_l_repressive_labor_pc L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc, fe vce(cluster cown)
est store lrag

gen lrag_sample = 0
replace lrag_sample = 1 if e(sample)

xtreg z2_polity_s_2010 L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc if lrag_sample==1, fe vce(cluster cown)
est store polity_lrsample

xtreg z2_polity_s_2010 z2_l_famfarms L4.z2_l_ln_tractors_ha ///
z2_l_log_gdp_per_cap_haber_men_2 z2_l_region_dems_pc if lrag_sample==1, fe vce(cluster cown)
est store famfarms_lrsample

***********
* TABLE 3 *
***********
esttab polity_cubic_time polyarchy_cubic_time libdem_cubic_time ///
polity_lrsample lrag famfarms_lrsample, ///
keep(L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc ///
z2_l_famfarms ///
z2_l_log_gdp_per_cap_haber_men_2 ///
z2_l_region_dems_pc) ///
order(L4.z2_l_ln_tractors_ha ///
z2_l_repressive_labor_pc ///
z2_l_famfarms ///
z2_l_log_gdp_per_cap_haber_men_2 ///
z2_l_region_dems_pc) ///
coeflabels(L4.z2_l_ln_tractors_ha "Tractors" ///
z2_l_repressive_labor_pc "LD Ag" ///
z2_l_famfarms "Fam Farms" ///
z2_l_log_gdp_per_cap_haber_men_2 "GDP/Cap" ///
z2_l_region_dems_pc "Reg Dem (\%)") ///
b(2) scalar(N_clust) aic se nogaps compress replace ///
star(* 0.10 ** 0.05 *** 0.01) ///
mtitles("Polity" "Polyarchy" "Lib Dem" "Polity" "Polity" "Polity") ///
addnotes("All models with country fixed effects, clustered standard errors" ///
"Models 3.1-3.3 with cubic time trend" ///
"Models 3.4-3.6 with Model 1.6 sample")
